﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Source
{
    public class FalseNearest
    {
        const long BOX = 1024;
        long[,] box;
        long length = long.MaxValue;
        long theiler = 0;
        int delay = 1;
        int maxdim = 5;
        int minemb = 1;
        int comp = 1;
        int maxemb = 5;
        long ibox = BOX - 1;
        long[] list;
        long toolarge;
        long[] vcomp;
        long[] vemb;
        double eps0 = 1.0e-5;
        double[][] series;
        double aveps, vareps;
        double varianz;

        private void mmb(long hdim, long hemb, double eps)
        {
            long i;
            long x, y;

            for (x = 0; x < BOX; x++)
                for (y = 0; y < BOX; y++)
                    box[x, y] = -1;

            for (i = 0; i < length - (maxemb + 1) * delay; i++)
            {
                x = (long)(series[0][i] / eps) & ibox;
                y = (long)(series[hdim][i + hemb] / eps) & ibox;
                list[i] = box[x, y];
                box[x, y] = i;
            }
        }

        byte findNearest(long n, int dim, double eps, double rf)
        {
            long x, y, x1, x2, y1, i, i1, ic, ie;
            long element, which = -1;
            double dx, maxdx, mindx = 1.1, hfactor, factor;

            ic = vcomp[dim];
            ie = vemb[dim];
            x = (long)(series[0][n] / eps) & ibox;
            y = (long)(series[ic][n + ie] / eps) & ibox;

            for (x1 = x - 1; x1 <= x + 1; x1++)
            {
                x2 = x1 & ibox;
                for (y1 = y - 1; y1 <= y + 1; y1++)
                {
                    element = box[x2, y1 & ibox];
                    while (element != -1)
                    {
                        if (Math.Abs(element - n) > theiler)
                        {
                            maxdx = Math.Abs(series[0][n] - series[0][element]);
                            for (i = 1; i <= dim; i++)
                            {
                                ic = vcomp[i];
                                i1 = vemb[i];
                                dx = Math.Abs(series[ic][n + i1] - series[ic][element + i1]);
                                if (dx > maxdx)
                                    maxdx = dx;
                            }
                            if ((maxdx < mindx) && (maxdx > 0.0))
                            {
                                which = element;
                                mindx = maxdx;
                            }
                        }
                        element = list[element];
                    }
                }
            }

            if ((which != -1) && (mindx <= eps) && (mindx <= varianz / rf))
            {
                aveps += mindx;
                vareps += mindx * mindx;
                factor = 0.0;
                for (i = 1; i <= comp; i++)
                {
                    ic = vcomp[dim + i];
                    ie = vemb[dim + i];
                    hfactor = Math.Abs(series[ic][n + ie] - series[ic][which + ie]) / mindx;
                    if (hfactor > factor)
                        factor = hfactor;
                }
                if (factor > rf)
                    toolarge++;
                return 1;
            }

            return 0;
        }

        public double[,] Calculate(string fileName, long dialogLenght, long exclude, int[] columns, long minEmbed, long maxEmbed, long delay, double rf, long theilerWindow, string separator)
        {
            double min, inter = 0.0, ind_inter, epsilon, av, ind_var;
            byte[] nearest;
            byte alldone;
            long i;
            int dim, emb;
            long donesofar;

            series = Routines.Routines.GetMultiSeries(fileName, out length, exclude, columns, separator);

            if (dialogLenght != long.MaxValue)
            {
                length = dialogLenght;
            }

            if (comp < columns.Length)
            {
                comp = columns.Length;
            }

            for (i = 0; i < comp; i++)
            {
                Routines.Routines.RescaleData(series[i], length, out min, out ind_inter);

                Routines.Routines.Variance(series[i], length, out av, out ind_var);

                if (i == 0)
                {
                    varianz = ind_var;
                    inter = ind_inter;
                }
                else
                {
                    varianz = (varianz > ind_var) ? ind_var : varianz;
                    inter = (inter < ind_inter) ? ind_inter : inter;
                }
            }

            list = new long[length];
            nearest = new byte[length];
            box = new long[BOX, BOX];

            vcomp = new long[maxdim];
            vemb = new long[maxdim];

            for (i = 0; i < maxdim; i++)
            {
                if (comp == 1)
                {
                    vcomp[i] = 0;
                    vemb[i] = i;
                }
                else
                {
                    vcomp[i] = i % comp;
                    vemb[i] = (i / comp) * delay;
                }

            }

            eps0 = 1.00 / 1023.00;

            double[,] result = new double[maxemb - minemb + 1, 4];

            for (emb = minemb; emb <= maxemb; emb++)
            {
                dim = emb * comp - 1;
                epsilon = eps0;
                toolarge = 0;
                alldone = 0;
                donesofar = 0;
                aveps = 0.0;
                vareps = 0.0;
                
                for (i = 0; i < length; i++)
                    nearest[i] = 0;

                while (alldone == 0 && (epsilon < 2.0 * varianz / rf))
                {
                    alldone = 1;

                    mmb(vcomp[dim], vemb[dim], epsilon);

                    for (i = 0; i < length - maxemb * delay; i++)
                        if (nearest[i] == 0)
                        {
                            nearest[i] = findNearest(i, dim, epsilon, rf);
                            
                            alldone &= nearest[i];

                            donesofar += (long)nearest[i];
                        }

                    epsilon *= Math.Sqrt(2.0);
                    
                    if (donesofar != 0)
                        eps0 = epsilon;
                }
                if (donesofar == 0)
                {
                    throw new Exception("Not enough points found!");
                }
                aveps *= (1.0 / (double)donesofar);
                vareps *= (1.0 / (double)donesofar);

                result[emb - minemb, 0] = dim + 1;
                result[emb - minemb, 1] = (double)toolarge / (double)donesofar;
                result[emb - minemb, 2] = aveps * inter;
                result[emb - minemb, 3] = Math.Sqrt(vareps) * inter;
            }

            return result;
        }
    }
}
